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ANALYSIS AND CONVERGENCE OF THE MAC SCHEME. 

1. THE LINEAR PROBLEM 


R. A. Nicolaides 1 
Department of Mathematics 
Carnegie Mellon University 
Pittsburgh, PA 15213 


ABSTRACT 

The MAC discretization of fluid flow is analyzed for the stationary Stokes equations. It 
is proved that the discrete approximations do in fact converge to the exact solutions of the 
flow equations. Estimates using mesh dependent norms analogous to the standard H and 
L 2 norms are given for the velocity and pressure respectively. 
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1. Introduction 


The original MAC (Marker and Cell) scheme of Harlow and Welch is now over thirty 
years old. Many evolutionary changes have occurred and the MAC approach is still a widely 
used finite difference method for incompressible flow computation. 

The MAC scheme is described in [Hirt 1979] and in most CFD books, for example 
[Fletcher 1989] and [Peyret and Taylor 1983]. It also forms the basis of several flow packages 
including SOLA [Hirt 1979] and SIMPLE [Patankar 1980]. 

More recently, it is finite element techniques which have dominated incompressible CFD. 
Although highly successful in many ways, the absence of fully satisfactory lower order 
schemes continues to be a stimulus to the development of new approaches. Recent work 
on finite elements may be found in [Girault and Raviart 1986], [Gunzburger 1989] and 
[Pironneau 1989]. 

The stability problems of low order finite elements [e.g Boland and Nicolaides 1985] 
suggest looking for stable schemes outside the finite element framework. The MAC scheme 
is one possibility. The standard MAC method suffers from being limited to rectangular 
meshes. A second difficulty is the apparent absence of any rigorous analysis for the stationary 
problem - the subject of this report. ([Porsching 1979] discusses the evolutionary problem.) 

On the first point, a generalization of the MAC ideas to triangular and other more 
general meshes was presented in [Nicolaides 1989]. In this complementary volume or covolume 
framework the MAC ideas are given a control volume interpretation. Complementary pairs 
of control volumes are used, and they can conveniently be taken as the triangles and polygons 
of a Delaunay- Voronoi mesh system. When specialized to a standard triangulation of, for 
example, a rectangular domain the precise MAC equations are recovered. [Nicolaides and 
Wu 1991] contains an implementation for a specific problem. [Hall, Cavendish, and Frey 
1991] contains a more general implementation. 

The covolume formulation is based on discretization of the vector operators grad, div, 
and curl and unlike the usual MAC derivation is coordinate free. The error estimates that 
are proved here follow from this formulation of the MAC scheme. A related approach was 
used in [Nicolaides 1991] to discretize and estimate some div-curl systems and several results 
from there will be used below. 

Before proceeding, a few technical comments are in order. First, we will present the 
analysis for two dimensional problems although the methods used are not inherently two 
dimensional. Second, the estimates for the velocity are in a discrete H 1 norm on Cl. This 
implies in particular, estimates for the vorticity on the boundary. While such estimates are 
desirable from a practical point of view, we have assumed that u> £ H 2 and p £ H 2 . This is 
more regularity than one would like to see. In our analysis, the MAC scheme itself is obtained 
as a discretization of the “strong” second order form of the equations. This already suggests 
the need for extra regularity. The third point is that while we use constant mesh spacings 
in the x and y directions, the domain need not be rectangular or convex. Coordinate free 
methods make this possible. Summation by parts formulas for example become inherently 
multidimensional and not limited to summations along grid lines. 
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2. Preliminaries 


For clarity, we assume a rectangular domain 0 and a primal Cartesian mesh having x 
and y spacings equal to h and h'. To avoid unnecessary complications, we will suppose 
that h = h'. It is easy to modify the results to cover the contrary case. In addition, by 
joining the centers of adjacent mesh rectangles (cells) a dual (“staggered”) mesh is made. 
At boundaries, nodes of the dual mesh are joined to the midpoints of the adjacent boundary 
edges. 

The N nodes of the primal mesh are numbered from left to right, and from bottom to 
top. The T nodes of the dual mesh are numbered in a similar way. There are T mesh cells in 
all, and there axe N dual mesh cells including the half size cells at the boundary. Similarly, 
the E edges of the primal mesh are labelled in some convenient way. There is a bijective 
correspondence between dual and primal edges, each edge crossing exactly one edge of the 
other set. The numbering of the dual edges reflects this, each being numbered the same as its 
primal companion. The cells, edges and nodes of the primal mesh are denoted generically by 
Tj, Gj and Vk respectively. Those of the dual mesh axe similarly denoted by primed quantities 
such as a'j. The lengths of the dual mesh edges are equal to h except near boundaries where 
they may be h/2. A direction is assigned to each primal edge according to the rule that 
positive is from low to high node number. The dual edges are directed by the convention 
that (<r£, a*) are oriented like the (x, y ) axes of the coordinate system. 

Defined on each primal edge aj is a component of the velocity field, directed in the positive 
direction of cr'. This velocity component is computed as n-u or its average, where u denotes 
the velocity vector at the midpoint of the primal edge, and n is the normal to the edge. 
Components of other vector fields are defined similarly. Such sets of normal components 
defined on the primal edges can be identified with R E . We will introduce an inner product 
into R E by 

{u,v)w := £ u j v jWj. 

<rj€ fi 

In this, Wj equals twice the area of the figure obtained by joining the nodes of a primal edge 
to the adjacent dual nodes. The sum is to be taken over all edges of the primal mesh. The 
associated norm is denoted by || • \\w- Clearly, it is twice a discrete L 2 norm. This inner 
product space is denoted by U, or by U 0 if the normal components assigned to the boundary 
edges are all zero. 

Various scalar fields including pressures are defined at the dual nodes. They can be 
identified with elements of R T . An inner product on R T is defined by 

(<f>, 0 )a := £ MA- 

TjCfi 

Here, the sum is over all mesh cells and A{ denotes the area of the i th cell. The associated 
norm is denoted by || • \\a- This inner product space is denoted by V. 

Analogous to V, scalar fields defined on the primal nodes can be identified with elements 
of R N , and an inner product defined by 

(V>> x)a' := £ TpkXkA' k , 

r'en 
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where the sum is over all dual cells including dual boundary cells, and A k denotes the area 
of the k th dual cell. The norm is denoted by || * \\ A * } and the inner product space by 5 , or 
by S 0 if the boundary values are all zero. Sometimes, it will be convenient to extend this 
norm and inner product over the interior or boundary nodes separately. A subscript will be 
used to denote such dependence thus: || • \\a\o- 

For each primal cell t*, discrete flux and divergence operators are defined on U by 

(Du)i:= £ u ih 

<rj£dri 


and 

(Du)i = (Du)i/Ai. 

By hj we mean hj negatively signed if the corresponding velocity component is directed 
towards the inside of r^, and positively signed otherwise. 

For each interior dual cell rj discrete circulation and curl operators are defined by 

( Cu)j := ^ u k h' k 


and 

(Cu)i = (<H /4 

This time the tilde produces a negative sign if the corresponding dual edge is directed against 
the positive sense of description of and a positive sign otherwise. 

This definition cannot be used in boundary dual cells. To extend it, we must require 
that tangential velocity components are specified along boundary segments defined by the 
intersections of consecutive dual mesh edges with T, Then we can define discrete circulations 
and curls in the same way even for the boundary dual cells. These extensions of C and C 
to the boundary are denoted by C \ and We will consider that the components of Cu 
are associated with interior nodes Vi and that the components of C\)U are associated with 
interior or boundary nodes as appropriate. 

Normal boundary components of u G U are denoted by tx|r- Sometimes, it will be 
convenient to use the same notation to include tangential components in the sense of the 
previous paragraph as well as normal components. It will be explicitly stated whenever such 
tangential boundary components are included. 

Slope operators are defined on V and S by 

:= 

where i\ and 12 are the nodes defining <j[ and the positive direction is from i\ to f 2 > an d 

f T> I \ ^ - V’il .• ^ 

{Rv) i •= — ^ — ) *2 > »i 

where <Tj is a primal edge with endpoints ii and ij. 

We will use the summation by parts formulas which follow: 

(t) (u, Rtp)w = (Cfeit, Vtt € U € S 
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where tangential components of u are zero. This notation is consistent with the definition of 
Cb only if the tangential components of u are zero. We will require a variant of this to cover 
the case when Rip is restricted to interior edges - those with at most one node on T. This is 

(i)' ( u , Rip)w = ( C b u , tp) A i Vu G U 0 'iip E S 

so that the normal boundary values of u are zero as well. 

(it) (u, G(p)w — -(Du, 4 >) a Vtz e Uq V<f> G V. 

These are easily proved by direct computation. 


3. Div-curl Results 


A covolume method and analysis was given in [Nicolaides 1991] for the div-curl system 

div u = p 
curl u = u) 

u • n| r - /. 

We will quote some results of [Nicolaides 1991] which are used below. 

In the notation introduced in section 2 the discrete div-curl approximation is 

Du = p 
Cu = u 
u|r = / 


where in each case the data is computed by simple averaging over the appropriate geometrical 
element. For example, 



and 



where <jj is numbered anticlockwise around T. 

It is worth remembering that Cu are the discrete curls (normalized circulations) taken 
around interior nodes only. 

If fl is multiply connected there is an additional condition [Nicolaides 1991], but we are 
considering rectangular domains here. The results stated below are valid generally. 

Define 


u 


( 2 ) 


■=T~nJ “• 

cr L Jo r' 


ids Ck 6 ft 


(i) 


where the integration is along the positive direction t of a' k . The superscript is chosen for 
consistency with [Nicolaides 1991]. If a k € T then we define 
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An error estimate follows: 


Theorem 1. Assume that the div-curl equations have a unique solution u with u G H 2 (Q). 
Then (a) the discrete div-curl equations have unique solution u, and (b) the estimate 

||u - uW\\ w < Kh 2 |u|h3(o) 

holds where K is independent of u and h. 

Proof. This is proved in the remarks following Theorem 6.1 in [Nicolaides 1991]. 

We will also need the following bound on the discrete solution by its data. 


Theorem 2. The solution of the system 


Du 

Cu 

u| r 


0 

u> 

0 


satisfies 

IM|w < K\\u\\ a , 

where K is a constant independent of ui and h. 

The norm on the right here refers to the summation over interior nodes of the primal 
mesh. 


Proof. See section 7 of [Nicolaides 1991]. 

4. Stokes Equations 

The inhomogeneous stationary Stokes equations are 

Au — VP = f in ft 

div u = 0 in ft 

u|r = g- 

Here, ft denotes a polygonal domain in R 2 , and T denotes its boundary. 

The standard weak form of the Stokes equations has unique solution u G H 1 (ft), p € 
Lg(ft) if f € L 2 (ft) and g G H 1 / 2 (T) [Girault and Raviart 1986]. 

The vorticity u) is defined to be curlu := d x v — d v u where the velocity field u := (u, u). 
In terms of u and using the incompressibility constraint, the momentum equation becomes 

—curl cu — VP = f (2) 
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where curl denotes the operator ( d v , -d x ). 

Let Oi denote an interior mesh segment, and let t and n be the unit tangent and unit 
normal directed along the positive directions of cr* and respectively. cr x may have at most 
one node on T. Let A and B denote s nodes where the positive direction is from A to B. 
Taking the inner product of the last equation with n and integrating along <j, from A to B 
gives 

— ^(B) + w(A) — J -^dt + J i-nds. ( 3 ) 

If, say, B is a boundary node we will use this formula to extend the vorticity to B. The 
extended value is well defined and independent of the particular line segment through B 
(or the point A on it) under the assumptions u; E VP € f E To 

prove this, denote the extension to B just defined by u>i(B), and let another one based on a 
positively directed line segment from C to B be called v 2 (B). Then it follows that 

MB) - uxfB)] + K A) - w(C) ] = f A + t ■„)*-/ & + { . n) * 

Jab an Jcb on 

so that 

[^{B) — u>i{B)\ — f (vp — f)-nds 

Jd{ABC) 1 

= 0 

where now n denotes the outer unit normal to the boundary d(ABC) of the triangle ABC 
and the trace theorem was used in association with (2). 

Note that since u> E by Sobolev’s lemma uj is uniformly continuous on Q and so 

can be uniquely extended to be continuous on fi. The extension of u> onto ft using ( 3 ) and 
the Sobolev extension coincide: for denoting by u>s(B) the Sobolev value it follows that with 
A E ft 

\u(B)-u s (B) | < |a;(J5)-a;(^)| + | Ws (B)- w (i)| 

< J AB \^\ds + \u s (B)-w s (A)\ 

< |v45| 1/2 ||p|| H j( n ) + \w s (B) — o> 5 (>1)| 

where the integration is performed along the line segment joining A and B , and Cauchy’s 
inequality and a trace theorem were used. The result follows from this. 


5. Discrete Stokes Equations 

The discretization of the Stokes equations is based on (3) obtained above. The momentum 
equation may be written as 


where 




,dp 1 f dp , 
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and n denotes the unit normal to <jj, and F denotes the vector of average normal components 
of f . In terms of the operators of section 2, the discrete equations are by definition 

-RJ -Gp' = F (?i e ft 
Du f = 0 
C h v! = a/. 

Boundary values for u are defined as simple averages of the normal and tangential com- 
ponents of g along the boundary edges of the mesh, so that a typical normal boundary value 

i s ^ 

v! = y / g * n da a » 6 T 

where n denotes the unit normal directed positively for the edge and h denotes the edge 
length. The typical tangential boundary value, associated to a primal boundary node v% is 
similarly the mean of the tangential component of g between the dual edge intersections 
enclosing 

To normalize the pressure approximation, we may impose 

(i,p')^ = o. 

The discrete momentum equation has one component per interior edge, and there is 
one incompressibility constraint for each primal cell. Since there is one unknown velocity 
component per interior edge and one pressure variable per cell there is a match between 
equations and unknowns, excluding the pressure normalization. 

It can be proved by direct calculation [Nicolaides 1988] that this approximation to the 
Stokes operator is precisely the MAC approximation. The equations would be identical were 
it not for our treatment of the data by exact integration in place of the quadrature likely to 
be used in reality. 

Explicit in the equations above is an unexpected connection with the velocity-vorticity 
formulation of the Stokes equations. From our discrete equations above we see that the 
MAC scheme can reasonably be described as a discretization of that formulation. But in the 
standard finite difference way, they are derived as an approximation to the primitive variable 
equations. So the MAC scheme simultaneously discretizes the two forms of the governing 
equations. Whatever advantages or disadvantages are perceived in either formulation are 
therefore present in the discrete scheme. From a slightly different viewpoint it follows that 
there are discrete analogs in the MAC framework (and in the covolume framework in general) 
of the transformations which enable us to go from one formulation to the other [Choudhury 
and Nicolaides 1991]. 

Including the pressure normalization there is one more equation than unknowns. This is 
not always convenient. We can avoid it by subtracting the mean pressure in the momentum 
equation. For this, let e € R T be the vector with ones in all its positions, and let A G R T 
have Aj := Aj/|J2|. In place of Gp we write G(I — eA l )p. The discrete Stokes equations 
become 

-RC b u' -G(I -eAy = F 
Du' = 0 
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together with the boundary equations. Now the pressure normalization is built in and the 
number of equations and unknowns and equations is the same. 

Theorem 1. The equations 


-RC b u' - G{I - eiV = F 

Du' = 0 

with ‘prescribed, normal and tangential boundary values it|r have a unique solution. 

Proof. Consider the homogeneous equations 

-RC b u' - G(I - eA l )p' = 0 

Du' — 0 

u'|r — 0. 

Included in the homogeneous boundary condition are both the normal and tangential values. 
Taking the inner product with u' € Uq we obtain 

(u 1 , RC b u')w + (u 1 , G(I — eA^p^w = 0. 

Using the summation formulas (i)' and ( ii ) reduces this to 

\\C b u'\\% = 0 

and in particular Cu' = 0. By Theorem 3.1 (a) it follows that u 1 = 0. Then G(I~eA t )p' = 0 
and (/ — eA l )p' is constant since the mesh is connected. But p' can differ from its mean 
by a constant only if the constant is zero, and so p' = 0. Uniqueness follows, and since the 
coefficient matrix is square so does existence. 


6. Error Bounds 
We will define 

(£>«)» := J ; u> dx dy i \ 6 

By Stokes’ theorem we also have 

(w„)i = -jj [ u ■ tds Vi 6 ft. 

S\ x Jdr{ 

Associated with Q v is a unique discrete velocity field v defined by 

Cv = Q v 
Dv = 0 
u|r = u' 


(4) 

(5) 
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where we recall that is defined as the simple average of the normal component g • n on 
each boundary edge. 

The field v defined in this way has no tangential components on T associated with it. 
We will now define them to be equal to those of u' . It is convenient to continue to call the 
augmented field v. Recall that the tangential components of u' on T are defined as averages 
of the prescribed tangential velocity. It follows that the tangential components of it v are 
all zero. 

The last equations are the standard covolume approximation to the system 

curl u = h) 
div u = 0 
u • n| r = g • n. 

The difference u' — v satisfies the discrete div-curl equations 

D( u' — v) = 0 
C(u' - v) = u/ — u> 

{y! — u)|r = 0 

and by Theorem 3.2, it follows that 

||xt / — w||w < -Kilo/ — WulU'- (6) 

In addition to the interior values, boundary values for u>„ eure defined by 

(&„)< := (C^); ViET. 

Subtracting the exact and approximate momentum equations and introducing p E V 
equal to the mean pressure in the primal mesh cells we obtain 

- R{cj' - Co v ) - G{p' -p) = —R(u> - a>„) - (m(^) - Gp ) cr E (7) 

The meaning of the notation a E ft is that the equation has one component for each edge of 
the primal mesh which is not on the boundary T. 

It follows that 

-(v! - v, R(u>' - a>„))w - (ix' - v, G(p' - p))w — 

dp 

(u 1 - v, R(oj - U>„))iy - (u 1 - V, — Gp)w £T 6 n 

and using both of the boundary conditions on v! — v and the summation formulas (i)' and 
(it) that 

—(Cb(u' — u), u' — & v )a' + (D(u' — u), p' — P)a — 

dp 

-(C b (u' - v), u - u v )a' - {u 1 - V, p(~) - Gp)w- 
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( 8 ) 


Thus, 

(u/ - u>„, a/ - Q v ) a , = (oj - Q v , uj - & V ) AI + (y! — v, /x(^) ~ Gp)w 
from which it follows by Cauchy’s inequality and bounding ||tx / — u||vv as in (6) that 

||a» — 6)4,11,4/ < ||u; — 6) v ||x' + — Gp^w 

In this equation, the A! norms are extended over all nodes including those on F. 

We will define a new field u* £ U as follows: on interior dual edges it is given by (1), and 
on dual boundary edges by a similar average. It also has tangential boundary values equal 
to those of u' . By definition of C{, it follows from (4) and (5) extended to T that 



u) gives the mean vorticity in each dual cell, including those on T. Notice that for interior 
dual cells u> v and u> coincide. 

Using (8) we now have with Ke{p) := ||/x(|E) - Gp\\ w 

ll"* — ©IU' < ||w' — uf v \\ A i + ||6>„ 

— \\u> — u) v \\ A ' -f \\u) — w v \\ A i + Ke(p) 

< \\v — G)\\ a < + 2\\u> — u) v \\ A i + Ke(p) 

< ||u> — u)\\ A ‘ + 2||u> — u>„||x',r + Kt{p) 

— ||u> — wU^'.n + ||u; — 6)||x',r + 2||a> — u^U^p + Ke(p). (9) 

The terms bounding the error can be estimated by approximation theory and Theorem 

1. 

6. Estimates 

To begin with we will estimate ||o> — 6)||x',n in (9). For this, consider the linear functional 

/ 1/2 f l /2 

/ u((,T})d(dr) 

-1/2 J -1/2 

on H 2 {Q X ), where Q x denotes the square -1/2 <( < 1/2, -1/2 < rj < 1/2. It is clear that 
the functional is bounded and is zero on linear functions. It follows that 

IM W )I ^ K\u\ 

Changing the variables in the integrals to x := £h and y := rjh gives eventually 

!(<*> - Q)j\ < r'j € fl. 

Squaring and summing over all rj £ fi gives 

||o/ — < K h 2 \w |^3(n)- 
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For h sufficiently small (< Mj7i(n)|/Mi7 a (n)l) it follows that 

||w - a>|U',n < Kh\u\H*(n ) 
which is the form we will use below. 

For the corresponding boundary term we consider the cell —1/2 < £ < 1/2, 0 < r} < 1/2. 
Then denoting this by Q 2 , we have 

|a,(0 ) 0)- 1 ^l/ Qa U,^^| < MO, 0)-M0, J)| + M°» 4 ) - ]q7\ Jq 3 u dr} \' 

The second term can be estimated by the same method used above for Q\. For the first 
term, 

I /^(O, ri)dr)\ < JTK( 0 , 011(0,*) 

J 0 

< K(\u>\hhq 3 ) + Mh’(Qj)) 

using the trace theorem. Changing variables to x := £h and y := rjh gives 

-i 

I / 4 w»(0, y)dy | < K(\w\ H i {Q )) + fc|w|H*(Q aifc )). 

Jo 

Squaring, multiplying by /i 2 and summing over the boundary nodes and incorporating the 
second boundary term yields 

M _ Mx'.r - # 2 (^ 4 Mff a (n) + ^ z Mff l (n)) 

so that for h sufficiently small we have the estimate 

||w — ^|U',r < Kh\u>\ H i(ny 

We shall amalgamate the interior and boundary estimates to get 

||w-w||x',n < Kh\v\ H i.(a) 

< Kh\u\n^(n) 


for h sufficiently small. 

Next we will estimate ||a>„ - 6>||^, r . The contribution to this sum associated with a fixed 
boundary node i has the form |(a> v )« - (^)«| 2 T- Suppressing the dependence on i we have 

|w v — w| = \(u* - v)x~ - (u* - v) 2 - + {u* — v) 3 h\/— 

< |(u* — u)i — (u* — v) 2 \/h + |(u* — u) 3 | / — . 

In this, the subscripts refer to the two dual edges of length £ and the interior dual edge of 
length h which are used to compute the discrete curl around the boundary node. There is 
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no contribution from the boundary itself because the tangential component of u* — v is zero 
there. The term (u* — u)i, for example, is the difference between the average normal velocity 
component measured along the normal from the boundary to the dual node of the boundary 
cell in question and the average of the (prescribed) normal velocity component measured 
along the boundary edge of the same boundary cell, (u* — v) 2 denotes the similar difference 
of averages on the opposite side of the node i. The union of these two boundary cells is the 
cell Q 2i h which was encountered in the previous estimate. 

To estimate the first term, we note that the functional (u* — v) x — (it* — v) 2 is bounded 
on H 2 (Q 3 ) and it can be checked that it vanishes on linear vector fields. Using the familiar 
argument and a scale change we obtain the estimate 


IK - w )i - ( u * - *>) 2 | < Kh\n\H7 {Qi h) . 

Summing the corresponding terms over T gives 

53 K u — u )i — ( u — u) 2 | 2 < -^ 2 ^ 2 | u Ih»(o)- 

r 

For the remaining part, using Theorem 3.1 we obtain 

£|( u '-v)|’ = £ l(“‘ - »)I* J 

r r 

n 

< K 2 h-W\vi\H*w 

< K 2 h 2 \u\ H 7(ny 

In this, the sum on the left is over all dual edges parallel to T appearing in the definition of 
the boundary circulations (in fact, the ‘dual boundary’ of fi) and the right hand sum in the 
second line is over interior edges of the primal mesh. Assembling the last two estimates we 
obtain 

ll^u — k’lU'.r < 

Estimation of the pressure term follows similar lines. The basic functional for this case 

,3(?) := c* *> - c L vitin - r:j> ^ ■ 

l 3 (-) is bounded on H 2 (Q 3 ), where Q 3 denotes the rectangle with corners at (— l,±l/2) and 
(l,±l/2), and is zero for linear functions. It follows that 


I ^3 (p) | < A|p|if2(Q s ). 


Introducing the scale changes x := (h, y := r/h we obtain 


1 f dp(s> V ) d _ Ph ~ Pfci 

h J<r k dx h 


< K\p\H>( Qs (h) 


12 



where the integration is along the positive direction of <7* and Qs{h) denotes the scaled region 
with cr* separating the primal mesh squares Tjt a and r^. A similar result holds when <7* is 
horizontal. Squaring and summing, it follows that 

m|£) - Gp\\w < Kh.\r\w m . (10) 

These estimates give the first part of the main result: 

Theorem 1. Under the regularity assumptions u G H 2 (n), u> G H 2 (Q) and p G H 2 (Q), the 
estimate 

IK - u>\ \ a< < Kh(\p\ H i(n) + Mif»(n)) 
holds for all h sufficiently small. 

Recalling that ||u||^ (n) = ||div u||£ 2(n) + ||curlu||£ 1(n) it follows that the norm of the 
error is indeed a discrete (mesh dependent) Hj(fl) norm. 


8. Pressure Error 

We begin this section by recalling the following standard result: 


Lemma. The equation 


divv = / £ lo(fi) 


has a solution v G Hj(fl) satisfying 


|v||tfi(n) < 


Proofs may be found in [Girault and Raviart 1986], and in [Temam 1984]. 

We will apply this result with 

/ := P - p' £ L o(ty 

where the right side denotes the piecewise constant function with these values in each primal 
mesh square. Clearly 

llp-p'IU’(n) = IIp-p'IU 

so that 

IMIifi(n) < K\\P-ii'\\a- (11) 

Next, introduce G Lio defined on each edge of the mesh as the mean of the normal 
component of v on the associated primal edge, i.e. 

:= r / V • n ds cr k G ft. 

h J < 7 * 
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Use of the divergence theorem shows that 


DvW = p - p' . ( 12 ) 

In addition, we have 

ll v(1) llw ^ -^ll v ll/f l (n) < K\\p - p'\\ A , (13) 

Only the first inequality is new. It is a consequence of the fact that the functional 

By := / vi(0, y)dy 

J —1/2 

is bounded on H 1 (Q4), where denotes the square with corners (±1/2, 0) and (0, ±1/2). 
Then changing the scale to h in both coordinate directions and summing over all the dual 
edges gives the result. 

Similar to we will need v* which is defined on the dual edges, as the mean of the 
tangential component along the dual edge in question: 

vx ‘ := ik\i ' vnd3 <€n - 

n points along o' in this. We will define tangential components for v* on the boundary 
segments delineated by dual mesh lines. These tangential values are defined to be zero. 

Now we have 

(CbV*y = — f v • t dt = — [ curlv dx dy 
T J Jd -i Tj Jr} 

from which it follows that 


\\C b v*\\ 2 AI < ||curlv||£ 2(n) < ||v||^i (n) < K\\p-p'\\ 2 A . (14) 

We will need an estimate for ||C , 6 V ^ 1 ^IU'* To obtain it, we first note the estimate 

||u* — < ■K'h||v|| H i (n) . (15) 

The technique for proving this is given in [Nicolaides 1991, Section 6]. We will omit the 
details here. Then 

IICI^’IU < ||C7 fc (v* - z.(‘>)|U- + IlCiu-IU- 

< jh" - ® (1) ||w + ||C»V*IU' 

< K\\p — p'\\ A (16) 

by (15), (11) and (14) . 

Now taking the inner product of i/ 1 ) with the basic error equation gives 

(® (1) . R(w' - w)V + G(p - P ')V = - Op) w 
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and using the summation formula (ii) gives 

\\P ~ p'Wa < ll^ ( 1 ) IU'IK - + ll u(1) MI/*(^) “ G P\\w- 

Then using (16) and (13) we finally obtain 

Up - p'\\a < K( ||o/ - w|L',n + ||p(|^) - Gp\\w- 

Using Theorem 6.1 and the approximation theory estimates we now have the second part 
of the main result: 

Theorem 1. Under the regularity assumptions u G H 2 (fi), a; G H 2 (Cl) and p G H 2 (£l), the 
estimate 

\\p' - P\\a < ^(|p|ff a (n) + l u lfl- J (0)) 
holds for all h sufficiently small. 

Thus the MAC scheme is first order accurate for the vorticity and pressure. It is not 
known yet whether the velocity is of higher order accuracy than the vorticity. Computations 
suggest that it is one order greater. [Nicolaides and Wu 1991]. 
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